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ABSTRACT 


We  are  developing  LITHOl.O:  an  updated  crust  and  lithosphere  model  of  the  Earth.  The  overall  plan  is  to  take  the 
popular  CRUST2.0  model — a  global  model  of  crustal  structure  with  a  relatively  poor  representation  of  the 
uppermost  mantle — and  improve  its  nominal  resolution  to  1  degree  and  extend  the  model  to  include  lithospheric 
structure.  The  new  model,  LITHOl.O,  will  be  constrained  by  many  different  datasets  including  extremely  large  new 
datasets  of  relatively  short  period  group  velocity  data.  Other  data  sets  include  (but  are  not  limited  to)  compilations  of 
receiver  function  constraints  and  active  source  studies. 

To  date,  we  have  completed  the  compilation  of  extremely  large  global  datasets  of  group  velocity  for  Rayleigh  and 
Love  waves  from  lOmHz  to  40mHz  using  a  cluster  analysis  technique.  We  have  also  extended  the  method  to 
measure  phase  velocity  and  are  complementing  the  group  velocity  with  global  data  sets  of  longer  period  phase  data 
that  help  to  constrain  deep  lithosphere  properties. 

To  model  these  data,  we  require  a  starting  model  for  the  crust  at  a  nominal  resolution  of  1  degree.  This  has  been 
developed  by  constructing  a  map  of  crustal  thickness  using  data  from  receiver  function  and  active  source 
experiments  where  available,  and  by  using  CRUST2.0  where  other  constraints  are  not  available.  Particular  care  has 
been  taken  to  make  sure  that  the  locations  of  sharp  changes  in  crustal  thickness  are  accurately  represented.  This  map 
is  then  used  as  a  template  to  extend  CRUST2.0  to  1  degree  nominal  resolution  and  to  develop  starting  maps  of  all 
crustal  properties. 

We  are  currently  modeling  the  data  using  two  techniques.  The  first  is  a  linearized  inversion  about  the  3D  crustal 
starting  model.  Note  that  it  is  important  to  use  local  eigenfunctions  to  compute  Frechet  derivatives  due  to  the 
extreme  variations  in  crustal  structure.  Another  technique  uses  a  targeted  grid  search  method.  A  preliminary  model 
for  the  crustal  part  of  the  model  will  be  presented. 
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OBJECTIVES 


The  main  objective  is  to  take  the  popular  CRUST2.0  model  -  a  global  model  of  crustal  structure  with  a  relatively 
poor  representation  of  the  uppermost  mantle  -  and  improve  its  nominal  resolution  to  1  degree  and  extend  the  model 
to  include  lithospheric  structure.  The  new  model,  LITHOl.O,  will  be  constrained  by  many  different  datasets 
including  extremely  large  new  datasets  of  relatively  short  period  group  velocity  data  measured  with  a  new  technique 
which  are  sensitive  to  lid  properties  as  well  as  crustal  thickness  and  average  crustal  properties.  Other  data  sets 
include  (but  are  not  limited  to)  compilations  of  receiver  function  constraints  and  active  source  studies. 

RESEARCH  ACCOMPLISHED 


Global  Surface  Wave  Group  Velocity  and  Phase  Velocity  Data  Sets 

In  a  previous  paper  (Masters  et  al.,  2010),  we  described  a  new  technique  to  measure  the  relative  group  arrival  times 
of  surface  wave  packets  that  allows  very  efficient  construction  of  large  data  sets.  We  have  applied  the  technique  to 
all  LH  data  at  the  Incorporated  Research  Institutions  for  Seismology  Data  Management  Center  (IRIS  DMC)  from 
events  larger  than  Ms=5.5  between  1976  and  2007.  We  also  include  some  BH  data  from  Program  for  the  Array 
Seismic  Studies  of  the  Continental  Lithosphere  (PASSCAL)  experiments  in  the  southern  hemisphere  and  the 
Portable  Observatories  for  Lithospheric  Analysis  and  Research  Investigating  Seismicity  (POLARIS)  network  to 
increase  ray  coverage.  The  Rayleigh  wave  group  velocity  data  set  is  complete  from  lOmHz  to  40mHz.  Similar 
datasets  for  Love  waves  are  now  complete. 


Figure  1.  Screen  capture  showing  the  processed  records  after  clustering — the  cluster  tree  is  on  the  right.  A 
vertical  cursor  is  used  to  select  the  degree  of  clustering — in  this  case,  there  is  one  useful  cluster 
comprising  about  80%  of  all  records.  The  top  records  which  are  noisy  form  a  set  of  small  clusters 
which  are  discarded.  This  example  is  for  a  period  of  100  seconds. 

At  very  long  periods  (greater  than  100s),  group  velocity  measurements  are  less  reliable  than  at  shorter  periods  but  it 
is  possible  to  modify  the  group  velocity  measurement  algorithm  to  measure  phase  velocity.  When  measuring  phase, 
we  are  no  longer  working  with  envelopes  so  we  must  be  very  careful  to  avoid  cycle-skipping.  We  have  found  that 
the  following  algorithm  works.  As  with  the  group  velocity  measurements,  we  band-pass  the  records  around  a 
frequency  of  interest  but  now  we  also  “undisperse”  the  records  using  a  ID  phase  velocity  curve  so  that  all  records 
are  arriving  at  a  nominal  distance  of  90  degrees.  We  correct  for  source  phase  and  amplitude  and  we  also  correct  for 
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the  predicted  phase  shift  due  to  3D  structure  by  using  a  nominal  phase  velocity  map.  This  latter  step  is  important  at 
periods  shorter  than  1 00  seconds  to  avoid  cycle  skipping  but  is  not  necessary  at  longer  periods.  The  algorithm  now 
follows  the  group  velocity  algorithm  where  we  cross-correlate  each  record  to  every  other  for  an  event  then  use  the 
cross-correlation  coefficients  to  perform  a  cluster  analysis  (Figure  1). 

Once  clusters  of  records  are  identified,  we  redo  the  cross-correlation  analysis  on  each  cluster  to  obtain  an  optimal  set 
of  relative  arrival  times.  Note  that  is  is  important  to  limit  how  much  each  record  can  be  time-shifted  in  the  cross¬ 
correlation  analysis.  We  have  found  that  we  must  limit  time  shifts  to  less  than  about  half  a  cycle.  This  is  why  it  is 
important  to  pre-shift  records  using  a  preliminary  phase  velocity  map.  We  start  the  analysis  at  long  periods  (200 
seconds  or  5mHz)  then  work  up  in  ImHz  increments  to  6mHz,  7mHz,  etc  using  the  phase  velocity  map  from  the 
previous  frequency  as  a  starting  phase  velocity  map  for  the  next  frequency.  Note  that  we  only  have  to  go  through  the 
process  of  selecting  clusters  at  5mHz  intervals.  For  intermediate  frequencies,  we  can  use  the  same  clusters  as 
determined  for  the  lower  frequency  analysis.  Thus,  the  6mHz  analysis  is  done  completely  automatically  using  the 
clusters  determined  in  the  5mHz  analysis  We  have  analyzed  Rayleigh  waves  from  5  to  25mHz  and  see  no  difficulty 
in  going  to  higher  frequencies.  To  date,  we  have  analyzed  data  from  1988  to  2004  which  results  in  over  310,000 
measurements  at  lOmHz. 

The  relative  arrival  times  are  inverted  for  a  phase  velocity  map  which  is  then  used  to  convert  the  data  to  absolute 
times.  That  is,  we  solve  for  a  mean  phase  velocity  and  a  time  shift  for  each  cluster  after  correcting  the  relative  arrival 
time  data  for  the  predictions  of  the  preliminary  phase  velocity  map.  The  absolute  times  are  then  used  to  determine 
the  final  phase  velocity  map.  Figure  2  illustrates  the  phase  velocity  map  derived  at  a  period  of  50  seconds  from 
295,000  phase  measurements.  Such  maps  provide  an  excellent  fit  to  the  data  with  a  variance  reduction  of  95% — but, 
more  importantly,  the  SMAD  of  the  phase  arrival  time  residuals  is  only  3.4  seconds  which  is  comparable  to  our 
measurement  error  (Figure  3).  We  conclude  that  our  final  dataset  of  phase  measurements  will  likely  be  even  larger 
than  those  for  group  velocity  and  appear  to  be  of  higher  precision.  We  are  currently  in  the  process  of  performing  a 
comprehensive  error  and  resolution  analysis  for  the  Rayleigh  wave  datasets. 


Figure  2.  Phase  velocity  map  at  50  seconds,  derived  from  295,000  phase  measurements. 
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Figure  3.  Histograms  showing  the  fit  of  the  phase  velocity  map  at  20mHz  to  the  phase  arrival  time  residuals. 

There  is  a  95%  variance  reduction  after  inversion  (left)  and  the  residual  SMAD  after  fitting  is  about 
3.4  seconds. 


Figure  4.  Crustal  thickness  (including  ice  and  water  layers)  from  a  compilation  of  many  studies  (see  text). 


Constructing  a  Starting  Crustal  Model 

As  a  precursor  to  getting  the  final  crust  and  lithosphere  model,  we  need  a  new  initial  crustal  model  with  1  degree 
nominal  resolution  that  accurately  matches  continental  shelves,  etc.  This  new  model,  CRUST1.0,  differs  from 
CRUST2.0  (Bassin  et  al.,  2000)  in  that  we  use  the  concept  of  “crustal  type”  to  specify  only  the  properties  of  the 
crystalline  crust  (Mooney  et  al.,  1998).  Topography  and  bathymetry  are  taken  from  ETOPOl  (NOAA’s  NGDC: 
Amante  and  Eakins,  2009).  The  ETOPOl  package  also  includes  ice  thickness  maps  for  Antarctica  and  Greenland 
both  of  which  we  included.  Information  of  sediment  structure  comes  from  the  three-layer  sediment  model  of  Laske 
and  Masters  (1997)  though  with  modifications  to  some  sediment  velocities.  We  also  specify  crustal  thickness 
independently  of  crustal  type.  The  crustal  thickness  map  has  been  compiled  from  a  literature  search  of  active  source 
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and  receiver  function  studies,  and  by  compiling  maps  of  crustal  thickness  taken  from  the  literature  (Artemieva, 
2006,  Tesauro  et  ah,  2007,  Grad  et  ah,  2008). 

Where  no  other  information  is  available,  we  have  used  crustal  thickness  maps  based  on  gravity  (e.g.,  most  of 
Antarctica)  or  we  have  used  the  crustal  thickness  of  CRUST2.0.  The  result  of  this  is  shown  in  Figure  4.  We  have 
also  assembled  an  initial  estimate  of  the  error  in  the  crustal  thickness  map  -  typically,  normal  oceanic  crustal 
thickness  is  likely  to  vary  by  less  than  2km,  where  we  have  direct  observational  constraints  on  continental  crustal 
thickness  we  use  an  uncertainty  of  3km,  where  we  have  a  map  from  the  literature  we  assume  an  uncertainty  of  5km, 
and  everywhere  else  (including  maps  based  on  gravity  interpretations)  we  assume  an  uncertainty  of  10km. 


Figure  5.  The  top  shows  the  observed  map  of  group  velocity  perturbations  (in  percent)  at  a  frequency  of 
30mHz.  The  bottom  shows  the  prediction  of  the  new  crustal  modal. 
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The  crystalline  crust  is  represented  by  three  layers.  We  have  found  that  a  total  of  32  crustal  types  gives  a  reasonable 
representation  of  the  crystalline  crust  in  CRUST  1.0  though  comparisons  with  the  data  suggest  some  additional 
crustal  types  might  be  appropriate. 

It  is  interesting  to  compare  the  predictions  of  the  new  starting  crustal  model  with  some  of  the  group  velocity  data 
sets.  Figure  5  (top)  shows  the  observed  group  velocity  map  at  30mHz  while  the  bottom  shows  the  predictions  of  the 
new  crustal  model  embedded  on  top  of  AK135.  Note  that  a  separate  group  velocity  calculation  is  done  for  each  (1 
degree)  pixel  of  the  maps.  There  are  some  interesting  discrepancies,  e.g.,  the  observed  map  suggests  that  we  should 
include  a  new  crustal  type  for  young  ocean  and  back  arc  basins.  There  are  also  some  interesting  discrepancies  where 
there  are  intermediate  thickness  sediment  piles  such  as  off  the  east  coast  of  the  United  States  of  America  and  the 
west  coast  of  Australia.  However,  we  feel  that  the  starting  model  is  sufficiently  close  for  our  purposes.  This  is 
verified  by  comparing  spherical  harmonic  expansions  of  the  two  maps  (Figure  6). 


Figure  6.  The  top  shows  the  amplitude  spectra  of  the  group  velocity  map  (red)  and  the  predicted  map  (blue) 
as  a  function  of  harmonic  degree  up  to  degree  180  (appropriate  for  a  1  degree  pixel  resolution).  The 
bottom  shows  the  correlation  between  the  two  along  with  the  99%  and  90%  confidence  levels. 

Clearly,  the  amplitude  spectra  of  the  two  maps  are  very  similar,  which  suggests  that  the  choice  or  regularization 
when  making  the  group  velocity  map  is  appropriate.  Furthermore,  the  correlation  is  extraordinarily  high  up  to  the 
highest  harmonic  degrees. 

Inverting  for  New  Crust  and  Lithosphere  Structure 

We  plan  to  use  two  techniques  to  obtain  improved  crustal  models.  The  first  technique  is  computationally  efficient 
but  relies  on  the  fact  that  we  have  a  starting  model  that  is  sufficiently  close  to  the  final  model  that  an  iterative, 
linearized  inversion  will  work.  To  do  this,  we  need  to  compute  group  velocity  kernels  for  each  pixel.  We  have  begun 
experimentation  changing  the  mean  crustal  velocity,  the  moho  depth,  and  the  mean  lid  velocity,  keeping  all  other 
things  fixed.  The  second  technique  uses  a  targeted  grid  search  method  that  was  presented  at  the  SSA  meeting  this 
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year  (Pasyanos  et  al.,  201 1).  This  technique  does  not  perform  any  linearization  but  is  computationally  much  more 
demanding  than  the  linearization  approach. 


Using  this  approach,  we  use  our  input  datasets  of  the  CRUST  1.0  crustal  profile,  crustal  thickness,  sediment 
thickness,  and  the  upper  mantle  model  LLNL-G3D  (Myers  et  al.,  2011;  Simmons  et  al.,  201 1)  to  create  an  initial 
model.  We  then  perturb  a  number  of  model  parameters  (crustal  velocities,  mantle  velocities,  crustal  thickness,  lid 
thickness)  and  create  a  suite  of  models.  The  models  are  then  tested  to  find  the  model  that  best  fits  the  surface  wave 
data.  For  example,  Figure  7  shows  an  example  of  a  grid  search  for  a  node  in  Siberia.  The  gray  lines  show  the 
dispersion  predictions  from  the  suite  of  models,  where  red  lines  indicate  better  fitting  models.  The  green  line,  which 
produces  the  best  fit  to  the  data  indicated  by  the  green  symbols,  is  shown  in  the  depth  profile  to  the  right.  This 
particular  node  required  a  slower  mantle  with  a  thick  lid. 


rcKle=  13664 1^=69.8517  ton=73.9D56 


Figure  7.  Example  of  targeted  grid  search  method  for  a  node  in  Siberia.  The  left  panel  shows  dispersion  for 
the  suite  of  models,  while  the  right  panel  shows  initial  and  new  model. 
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Figure  8.  Crustal  thickness  maps  of  a  global,  low-resolution  model  and  a  regional,  high-resolution  model. 


125 

Approved  for  public  release;  distribution  is  unlimited. 


2011  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


We  are  now  in  a  position  to  assemble  the  results  from  individual  nodes  to  start  and  construct  the  LITHOl.O  model. 
Figure  8  shows  a  global,  low-resolution  model  where  the  inversion  were  performed  at  tessellated  nodes  with  about 
8°  spacing  alongside  a  high-resolution  regional  model  of  the  Middle  East  where  the  node  spacing  is  about  1°. 

CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  completed  the  construction  of  a  new  large  global  dataset  of  Rayleigh  and  Love  wave  group  arrival  times 
using  a  new  efficient  measurement  technique  and  have  evaluated  the  data  by  constructing  group  velocity  maps  using 
simple  ray  theory.  We  have  extended  the  technique  to  measure  phase  and  the  new  datasets  of  phase  measurements 
are  likely  to  be  as  large  and  more  precise  than  the  group  velocity  data  sets. 

We  have  constructed  a  new  crustal  model  with  1  degree  nominal  resolution  as  a  starting  point  for  inversion.  This 
model  predicts  group  velocity  perturbations  with  very  similar  spectral  characteristics  to  the  data  though  detailed 
comparisons  suggests  further  improvement  of  the  starting  model  by  adding  a  few  new  crustal  types  and  by  adjusting 
some  sediment  velocities. 

We  are  now  in  a  position  to  begin  detailed  modeling  of  the  data  sets  for  the  new  crust  and  lithosphere  model. 
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